************************************************
* Table A8: Multiple hypothesis test corrections
************************************************


***************************************
**** PRE-AMBLE ****
***************************************

clear all
set more off
pause on
// version 16 

qui adopath ++ "$dir/ado"

global WyoungBootstraps $reps
global seed 4832

***************************************
**** 	LOAD MHT FUNCTIONS 		****
***************************************

*Anderson Code converted into a function that reads in the p-values stored in
*matrix P
cap program drop qValues
program define qValues
  syntax, MATRIX(string)
  *Anderson Code 
* This code generates BKY (2006) sharpened two-stage q-values as described in Anderson (2008), "Multiple Inference and Gender Differences in the Effects of Early Intervention: A Reevaluation of the Abecedarian, Perry Preschool, and Early Training Projects", Journal of the American Statistical Association, 103(484), 1481-1495

* BKY (2006) sharpened two-stage q-values are introduced in Benjamini, Krieger, and Yekutieli (2006), "Adaptive Linear Step-up Procedures that Control the False Discovery Rate", Biometrika, 93(3), 491-507

* Last modified: M. Anderson, 11/20/07
* Test Platform: Stata/MP 10.0 for Macintosh (Intel 32-bit), Mac OS X 10.5.1
* Should be compatible with Stata 10 or greater on all platforms
* Likely compatible with with Stata 9 or earlier on all platforms (remove "version 10" line below)

version 10

****  INSTRUCTIONS:
****    Please start with a clear data set
****  When prompted, paste the vector of p-values you are testing into the "pval" variable
****  Please use the "do" button rather than the "run" button to run this file (if you use "run", you will miss the instructions at the prompts)

pause on
set more off
save "TEMP.dta", replace 
clear 

if _N>0 {
  display "Please clear data set before proceeding"
  display "After clearing, type 'q' to resume"
  pause
  } 

*quietly gen float pval = .
svmat `matrix'
rename `matrix'1 pval

display "***********************************"
display "Please paste the vector of p-values that you wish to test into the variable 'pval'"
display "After pasting, type 'q' to resume"
display "***********************************"

*pause

* Collect the total number of p-values tested

quietly sum pval
local totalpvals = r(N)

* Sort the p-values in ascending order and generate a variable that codes each p-value's rank

quietly gen int original_sorting_order = _n
quietly sort pval
quietly gen int rank = _n if pval~=.

* Set the initial counter to 1 

local qval = 1

* Generate the variable that will contain the BKY (2006) sharpened q-values

gen bky06_qval = 1 if pval~=.

* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.


while `qval' > 0 {
  * First Stage
  * Generate the adjusted first stage q level we are testing: q' = q/1+q
  local qval_adj = `qval'/(1+`qval')
  * Generate value q'*r/M
  gen fdr_temp1 = `qval_adj'*rank/`totalpvals'
  * Generate binary variable checking condition p(r) <= q'*r/M
  gen reject_temp1 = (fdr_temp1>=pval) if pval~=.
  * Generate variable containing p-value ranks for all p-values that meet above condition
  gen reject_rank1 = reject_temp1*rank
  * Record the rank of the largest p-value that meets above condition
  egen total_rejected1 = max(reject_rank1)

  * Second Stage
  * Generate the second stage q level that accounts for hypotheses rejected in first stage: q_2st = q'*(M/m0)
  local qval_2st = `qval_adj'*(`totalpvals'/(`totalpvals'-total_rejected1[1]))
  * Generate value q_2st*r/M
  gen fdr_temp2 = `qval_2st'*rank/`totalpvals'
  * Generate binary variable checking condition p(r) <= q_2st*r/M
  gen reject_temp2 = (fdr_temp2>=pval) if pval~=.
  * Generate variable containing p-value ranks for all p-values that meet above condition
  gen reject_rank2 = reject_temp2*rank
  * Record the rank of the largest p-value that meets above condition
  egen total_rejected2 = max(reject_rank2)

  * A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
  replace bky06_qval = `qval' if rank <= total_rejected2 & rank~=.
  * Reduce q by 0.001 and repeat loop
  drop fdr_temp* reject_temp* reject_rank* total_rejected*
  local qval = `qval' - .001
}
  

quietly sort original_sorting_order
pause off
*set more on

display "Code has completed."
display "Benjamini Krieger Yekutieli (2006) sharpened q-vals are in variable 'bky06_qval'"
display "Sorting order is the same as the original vector of p-values"

mkmat pval, matrix(PVAL)
mkmat bky06_qval, matrix(Q)
mat list PVAL 
mat list Q
use "TEMP.dta", clear
cap erase "TEMP.dta"

end

*Create multiple hypothesis testing program
cap program drop MHT
program define MHT
  syntax, MODELS(string asis) TESTS(varlist) BOOTSTRAPS(integer) STRATA(varname) FILENAME(string) PANEL(string)

  save "$dir/TEMP1.dta", replace 

  *Create locals and new variables to store output 
  tokenize `"`models'"'
  local nmodels = wordcount("`tests'")
  
  *Westfall-Young with strata
  wyoung, cmd(`models') familyp(`tests') bootstraps(`bootstraps') cluster(pid) seed($seed) strata(`strata') replace 
  drop psidak
  order pwyoung, last 
  save "$dir/TEMP2.dta", replace 

  /*Westfall-young without strata 
  use "$dir/TEMP1.dta", clear
  wyoung, cmd(`models') familyp(`tests') bootstraps(`bootstraps') cluster(pid) seed($seed) replace 
  keep k pwyoung
  rename pwyoung pwyoung_no_strata
  merge 1:1 k using "$dir/TEMP2.dta"
  drop _merge 
  order pwyoung pwyoung_no_strata, last 
  save "$dir/TEMP2.dta", replace 
  */

  *Anderson Q-value correction
  use "$dir/TEMP1.dta", clear
  matrix P = J(`nmodels', 1, .)
  forvalues i = 1/`nmodels'{
    di "`i'"
    local var = word("`tests'", `i')   
    ``i''
    test `var'
    matrix P[`i', 1] = `r(p)'
  }

  qValues, matrix(P)
  use "$dir/TEMP2.dta", replace 
  svmat Q 
  svmat PVAL
  assert abs(PVAL1 - p) < 0.00001  
  drop PVAL1
  rename Q1 fdr_qvalue

  *Output latex table 
  gen col1 = . 
  gen col2 = . 
  order col1 col2 coef-fdr_qvalue, last 
  mkmat col1 - fdr_qvalue, matrix(A)
  frmttable using "$oa/$t/Table_A8/`filename'", statmat(A) sdec(3) ///
  ctitle("\uline{\hfill \textbf{`panel'} \hfill}", "",  "", "", "", "", "", "", "", \, "Model" "Variable" "Coefficient" "SE" "p-value" "\thead{Bonferroni \\ p-value}" "\thead{Westfall-Young \\ p-value}" ///
         "\thead{FDR \\ q-value}" ) multicol(1, 1, 8) tex fragment nocenter  replace
  end 


**********************************************
****    Figure - Thoughts while working
**********************************************

use "$dir/data/analysis_main.dta", clear

#delimit ; 
MHT, 
  models( 
    "reghdfe thinktask cash loan_moneylender numloans_worry bs_labor_30_nm bs_labor_30_miss i.future_fin_worry_nm if pid1==1 & lateprime==0 , vce(robust) a(roundgroup)"
    "reghdfe thinkonlytask cash loan_moneylender numloans_worry bs_labor_30_nm bs_labor_30_miss i.future_fin_worry_nm if pid1==1 & lateprime==0 , vce(robust) a(roundgroup)")
tests(
    cash
    cash)
  bootstraps($WyoungBootstraps)
  strata(roundgroup)
  filename(fig3_pval_corrections)
  panel("PANEL A: Worries (Figure III)");
#delimit cr 


***************************************
****   EXPENDITURE  ****
***************************************


use "$dir/data/analysis_expenditure.dta", clear

foreach x in loancredit_any {
  gen  `x'_temp1 =`x' if postpay_3d==1
  egen `x'_temp2 = max(`x'_temp1), by(pid)
}

foreach x in loancredit tothhexpense foodexp dressexp dailyexp medicexp tobaccoexp totexpense_top{
  gen `x'_temp1 =`x' if postpay_3d==1
  egen `x'_temp2 = total(`x'_temp1), by(pid)
}

#delimit ; 
MHT, 
  models( 
    "areg loancredit_temp2 cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss if pid1==1, cluster(pid) a(roundgroup)"
    "areg loancredit_any_temp2 cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss if pid1==1, cluster(pid) a(roundgroup)"
    "areg tothhexpense_temp2 cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss if pid1==1, cluster(pid) a(roundgroup)"
    "areg foodexp_temp2 cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss if pid1==1, cluster(pid) a(roundgroup)"
    "areg dressexp_temp2 cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss if pid1==1, cluster(pid) a(roundgroup)"
    "areg dailyexp_temp2 cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss if pid1==1, cluster(pid) a(roundgroup)"
    "areg medicexp_temp2 cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss if pid1==1, cluster(pid) a(roundgroup)"
    "areg tobaccoexp_temp2 cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss if pid1==1, cluster(pid) a(roundgroup)"
    "areg totexpense_top_temp2 cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss if pid1==1, cluster(pid) a(roundgroup)")
tests(
    cash
    cash
    cash
    cash
    cash
    cash
    cash
    cash
    cash 
    )
  bootstraps($WyoungBootstraps)
  strata(roundgroup)
  filename(table2_pval_corrections)
  panel("PANEL B: Expenditure (Table II)");
#delimit cr 



***************************************
**** 		PRODUCTIVITY 		****
***************************************

use "$dir/data/analysis_main.dta", clear

* Heterogeneity continuous
cap drop *_hiw
foreach x of varlist cash_postpay cash_announce{
  gen `x'_wealth_ind_avg = `x'*wealth_ind_avg
} 

* Heterogeneity binary
cap drop *_hiw
foreach x of varlist cash_postpay cash_announce{
  gen `x'_wealth_hi_avg = `x'*wealth_hi_avg
}   

*PRODUCTION
#delimit ; 
MHT, 
  models( 
    "reghdfe n_productivity cash_postpay cash_announce blhrprod blhrprod_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay cash_announce blhrprod blhrprodsq blhrprod_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay cash_announce blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay cash_announce blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup day Hour_Start)"
    "reghdfe n_productivity cash_postpay cash_announce cash_postpay_wealth_ind_avg cash_announce_wealth_ind_avg wealth_ind_avg blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay cash_announce cash_postpay_wealth_ind_avg cash_announce_wealth_ind_avg wealth_ind_avg blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay cash_announce cash_postpay_wealth_hi_avg cash_announce_wealth_hi_avg wealth_hi_avg blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay cash_announce cash_postpay_wealth_hi_avg cash_announce_wealth_hi_avg wealth_hi_avg blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)")
  tests(
    cash_postpay
    cash_postpay
    cash_postpay
    cash_postpay
    cash_postpay
    cash_postpay_wealth_ind_avg
    cash_postpay
    cash_postpay_wealth_hi_avg)
  bootstraps($WyoungBootstraps)
  strata(roundgroup)
  filename(table3_pval_corrections)
  panel("PANEL A: Production (Table III)");
#delimit cr 





***************************************
****    ATTENTIVENESS     ****
***************************************


use "$dir/data/analysis_main.dta", clear

* Heterogeneity continuous
cap drop *_hiw
foreach x of varlist cash_postpay cash_announce{
  gen `x'_wealth_ind_avg = `x'*wealth_ind_avg
} 

* Heterogeneity binary
cap drop *_hiw
foreach x of varlist cash_postpay cash_announce{
  gen `x'_wealth_hi_avg = `x'*wealth_hi_avg
}   
  
*PRODUCTION
#delimit ; 
MHT, 
  models( 
    "reghdfe att_ind cash_postpay cash_announce blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe hiatt_ind cash_postpay cash_announce blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe att_ind cash_postpay cash_announce cash_postpay_wealth_ind_avg cash_announce_wealth_ind_avg wealth_ind_avg blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe att_ind cash_postpay cash_announce cash_postpay_wealth_ind_avg cash_announce_wealth_ind_avg wealth_ind_avg blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe att_ind cash_postpay cash_announce cash_postpay_wealth_hi_avg cash_announce_wealth_hi_avg wealth_hi_avg blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe att_ind cash_postpay cash_announce cash_postpay_wealth_hi_avg cash_announce_wealth_hi_avg wealth_hi_avg blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe hiatt_ind cash_postpay cash_announce cash_postpay_wealth_hi_avg cash_announce_wealth_hi_avg wealth_hi_avg blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe hiatt_ind cash_postpay cash_announce cash_postpay_wealth_hi_avg cash_announce_wealth_hi_avg wealth_hi_avg blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)")
tests(
    cash_postpay
    cash_postpay
    cash_postpay
    cash_postpay_wealth_ind_avg
    cash_postpay
    cash_postpay_wealth_hi_avg
    cash_postpay
    cash_postpay_wealth_hi_avg)
  bootstraps($WyoungBootstraps)
  strata(roundgroup)
  filename(table4_pval_corrections)
  panel("PANEL B: Attention (Table IV)");
#delimit cr 


**********************************************
************    PIECERATE     *************
**********************************************


use "$dir/data/analysis_main.dta", clear

	
keep if pid1==1
keep pid blhrprod blhrprodsq blhrprod_miss blatt blattsq blatt_miss ///
	educ educsq married widowed divorced exper expersq wealth_hi_avg wealth_ind_avg hasnocredits hasloan loan_moneylender anyloans_worry richhouse richhouse_miss canget1k canget1k_miss ownland ownland_miss age_nm age_miss read_odiya_nm read_odiya_miss anychildren_nm children_miss caslabor_nm caslabor_miss bs_labor_7_nm bs_labor_7_miss bs_labor_30_nm bs_labor_30_miss future_fin_worry_nm worry1_nm worry2_nm future_fin_worry_miss bl_loanamt_nm bl_loanamt_miss anyloans_dueyr0_nm numloans_dueyr0_nm numloans_dueyr0_miss

sort pid
tempfile mainround
save `mainround', replace

use "$dir/data/analysis_all.dta", clear
replace piecerate = 0 if day==3 & pid==438
drop if pid==431

foreach x in blhrprod blhrprodsq blhrprod_miss blatt blattsq blatt_miss educ educsq married widowed divorced exper expersq wealth_hi_avg wealth_ind_avg hasnocredits hasloan loan_moneylender anyloans_worry richhouse richhouse_miss canget1k canget1k_miss ownland ownland_miss age_nm age_miss read_odiya_nm read_odiya_miss anychildren_nm children_miss caslabor_nm caslabor_miss bs_labor_7_nm bs_labor_7_miss bs_labor_30_nm bs_labor_30_miss future_fin_worry_nm worry1_nm worry2_nm future_fin_worry_miss bl_loanamt2_nm bl_loanamt2_miss bl_loanamt_nm bl_loanamt_miss anyloans_dueyr0_nm numloans_dueyr0_nm numloans_dueyr0_miss{
	capture confirm variable `x'
	if _rc==0 drop `x'
}

keep if round > 15
sort pid, stable
merge m:1 pid using `mainround'
drop if _merge==2

cap drop pid1 pidday1
egen pid1=tag(pid) if attend==1 & day!=1
egen pidday1=tag(pid day) 


*** wage and piece rate
replace baserate=. if attend==0
gen wage = baserate + piecerate*day_output   // does not include bonuses and game payments
replace wage=. if attend==0
gen lwage = log(wage)

gen pr0=piecerate==0
gen pr2=piecerate==2
gen pr3=piecerate==3
gen pr4=piecerate==4

label var pr0 "Piece rate = Rs. 0"
label var pr2 "Piece rate = Rs. 2"
label var pr3 "Piece rate = Rs. 3"
label var pr4 "Piece rate = Rs. 4"

gen lpiecerate = log(piecerate)
label var lpiecerate "Log(piece rate)"
label var piecerate "Piece rate"



* One person dropped out at the beginning in the main rounds; no baseline info
gen baseline_miss = pid==455
foreach x of varlist blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss blatt blatt_miss read_odiya_nm read_odiya_miss divorced{
	replace `x' = 0 if baseline_miss==1
}

* one person left early on one day
gen attend2 = attend
replace attend2 = 1 if pid==17 & day==4 & slot==5


save "TEMP3.dta", replace 

* conditional on attendance

#delimit ; 
MHT, 
  models( 
    "reghdfe n_productivity piecerate blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe n_productivity lpiecerate blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe n_productivity pr3 pr4 blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe n_productivity pr3 pr4 blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe att_ind piecerate blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe att_ind lpiecerate blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe att_ind pr3 pr4 blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe att_ind pr3 pr4 blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe attend piecerate blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & pidday1==1, vce(cl pid) absorb(round ) "
    "reghdfe attend lpiecerate blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & pidday1==1, vce(cl pid) absorb(round ) "
    "reghdfe attend pr3 pr4 blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & pidday1==1, vce(cl pid) absorb(round ) "
    "reghdfe attend pr3 pr4 blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & pidday1==1, vce(cl pid) absorb(round )")
tests(
    piecerate
    lpiecerate
    pr3
    pr4
    piecerate
    lpiecerate
    pr3
    pr4
    piecerate
    lpiecerate
    pr3
    pr4)
  bootstraps($WyoungBootstraps)
  strata(round)
  filename(table5_pval_corrections)
  panel("PANEL C: Piece Rate (Table V)");
#delimit cr 


use "TEMP3.dta", clear 

#delimit ; 
MHT, 
  models( 
    "reghdfe n_productivity piecerate blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe n_productivity lpiecerate blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe n_productivity pr3 pr4 blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe n_productivity pr3 pr4 blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) ")
tests(
    piecerate
    lpiecerate
    pr3
    pr4)
  bootstraps($WyoungBootstraps)
  strata(round)
  filename(table5_pval_corrections_production)
  panel("PANEL A: Piece Rate Production (Table V)");
#delimit cr 

use "TEMP3.dta", clear 

#delimit ; 
MHT, 
  models( 
    "reghdfe att_ind piecerate blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe att_ind lpiecerate blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe att_ind pr3 pr4 blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) "
    "reghdfe att_ind pr3 pr4 blatt blatt_miss blhrprod_miss read_odiya_nm read_odiya_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & attend==1, vce(cl pid) absorb(round ) ")
tests(
    piecerate
    lpiecerate
    pr3
    pr4)
  bootstraps($WyoungBootstraps)
  strata(round)
  filename(table5_pval_corrections_attention)
  panel("PANEL B: Piece Rate Attention (Table V)");
#delimit cr 


use "TEMP3.dta", clear 

#delimit ; 
MHT, 
  models( 
    "reghdfe attend piecerate blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & pidday1==1, vce(cl pid) absorb(round ) "
    "reghdfe attend lpiecerate blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & pidday1==1, vce(cl pid) absorb(round ) "
    "reghdfe attend pr3 pr4 blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & pidday1==1, vce(cl pid) absorb(round ) "
    "reghdfe attend pr3 pr4 blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss baseline_miss if day>1 & piecerate>0 & pidday1==1, vce(cl pid) absorb(round )")
tests(
    piecerate
    lpiecerate
    pr3
    pr4)
  bootstraps($WyoungBootstraps)
  strata(round)
  filename(table5_pval_corrections_attendance)
  panel("PANEL C: Piece Rate Attendance (Table V)");
#delimit cr 


cap erase "TEMP3.dta"



**********************************************
************    FAIRNESS CONCERNS     *************
**********************************************

use "$dir/data/analysis_main.dta", clear

#delimit ; 
MHT, 
  models( 
    "reghdfe n_productivity post_announce_1_cash post_announce_2_cash post_announce_3p_cash post_announce_1 post_announce_2 post_announce_3p blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss if pd==1, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity post_announce_1_cash post_announce_2_cash post_announce_3p_cash post_announce_1 post_announce_2 post_announce_3p blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss if pd==1, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity post_announce_1_cash post_announce_2_cash post_announce_3p_cash post_announce_1 post_announce_2 post_announce_3p cash_postpay blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity post_announce_1_cash post_announce_2_cash post_announce_3p_cash post_announce_1 post_announce_2 post_announce_3p cash_postpay blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity post_announce_1_cash post_announce_2_cash post_announce_3p_cash post_announce_1 post_announce_2 post_announce_3p cash_postpay blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity payday_cash payday papc_cash2 payday_cash_wb payday_wb blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss if pd==1, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity payday_cash payday papc_cash2 payday_cash_wb payday_wb blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss if pd==1, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity payday_cash payday papc_cash2 payday_cash_wb payday_wb blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss if pd==1, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity payday_cash payday papc_cash2 payday_cash_wb payday_wb cash_postpay blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity payday_cash payday papc_cash2 payday_cash_wb payday_wb cash_postpay blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity payday_cash payday papc_cash2 payday_cash_wb payday_wb cash_postpay blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity payday_cash payday papc_cash2 payday_cash_wb payday_wb cash_postpay blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)")
  tests(
   post_announce_1_cash 
   post_announce_2_cash
   post_announce_1_cash 
   post_announce_2_cash
   cash_postpay
   papc_cash2
   payday_cash 
   payday_cash_wb
   papc_cash2
   payday_cash
   payday_cash_wb
   cash_postpay)
  bootstraps($WyoungBootstraps)
  strata(roundgroup)
  filename(table6_pval_corrections)
  panel("PANEL D: Announcement Effects (Table VI)");
#delimit cr 


***************************************
****    NUTRITION     ****
***************************************

use "$dir/data/analysis_main.dta", clear


*** Last 2 hours - no decline in effects by last slot of the day

* shorten names
gen l1fs  = last1fullslot 
gen l2fs = last12fullslot
gen halfl1 = ((slot>=5 & round<=2) | (slot>=6 & round==3) | (Hour_Start==1 & round>3)) if round<=15
gen halfl2 = ((slot>=5 & round<=2) | (slot>=6 & round==3) | (Hour_Start==12 & round>3) | (Hour_Start==1 & round>3)) if round<=15

cap drop *_hiw

foreach x of varlist l1fs l2fs halfl1 halfl2 slot{
	foreach y of varlist cash_postpay {
		gen `y'_`x' = `y'*`x'
		gen `y'_`x'_hiw = `y'_`x'*wealth_hi_avg
	}
	gen `x'_hiw = `x'*wealth_hi_avg
}

foreach x of varlist cash_postpay cash_announce{
	gen `x'_hiw = `x'*wealth_hi_avg
}

save "TEMP3.dta", replace 

#delimit ; 
MHT, 
  models( 
    "reghdfe n_productivity slot cash_postpay_slot cash_postpay cash_announce blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity slot cash_postpay_slot cash_postpay cash_announce blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity slot cash_postpay_slot cash_postpay cash_announce slot_hiw cash_postpay_slot_hiw cash_postpay_hiw cash_announce_hiw wealth_hi_avg blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity slot cash_postpay_slot cash_postpay cash_announce slot_hiw cash_postpay_slot_hiw cash_postpay_hiw cash_announce_hiw wealth_hi_avg blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity slot cash_postpay_slot cash_postpay cash_announce slot_hiw cash_postpay_slot_hiw cash_postpay_hiw cash_announce_hiw wealth_hi_avg blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity slot cash_postpay_slot cash_postpay cash_announce slot_hiw cash_postpay_slot_hiw cash_postpay_hiw cash_announce_hiw wealth_hi_avg blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay_halfl2 cash_postpay cash_announce blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay_halfl2 cash_postpay cash_announce blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay_halfl1 cash_postpay cash_announce blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)"
    "reghdfe n_productivity cash_postpay_halfl1 cash_postpay cash_announce blhrprod blhrprodsq blhrprod_miss bs_labor_30_miss future_fin_worry_miss i.pd if pd>0, vce(cl pid) absorb(roundgroup)")
tests(
    cash_postpay
    cash_postpay_slot
    cash_postpay
    cash_postpay_slot
    cash_postpay_hiw
    cash_postpay_slot_hiw
    cash_postpay
    cash_postpay_halfl2
    cash_postpay
    cash_postpay_halfl1)
  bootstraps($WyoungBootstraps)
  strata(roundgroup)
  filename(table7_pval_corrections_production)
  panel("PANEL F: Nutrition Channel on Production (Table VII)");
#delimit cr 

use "TEMP3.dta", clear 

keep if round>3 & pid1==1

#delimit ; 
MHT, 
  models( 
    "reghdfe bf2d_avg cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss, cluster(pid) a(roundgroup)"
    "reghdfe rice2d_avg cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss, cluster(pid) a(roundgroup)"
    "reghdfe rice2dq cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss, cluster(pid) a(roundgroup)"
    "reghdfe veg2d_avg cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss, cluster(pid) a(roundgroup)"
    "reghdfe other2d_avg cash blhrprod_miss canget1k_miss ownland_miss bs_labor_7_miss bs_labor_30_miss future_fin_worry_miss, cluster(pid) a(roundgroup)")
  tests(
    cash    
    cash
    cash
    cash
    cash)
  bootstraps($WyoungBootstraps)
  strata(roundgroup)
  filename(table7_pval_corrections_breakfast)
  panel("PANEL E: Nutrition Channel Breakfast Measures (Table VII)");
#delimit cr 


cap erase "TEMP3.dta" 
cap erase "$dir/TEMP2.dta" 
cap erase "$dir/TEMP1.dta" 


